# d维空间有心势中粒子轨道的稳定性

import sympy as sp

r = sp.symbols('r',real=True,positive=True)
theta,pr,ptheta = sp.symbols('theta,pr,ptheta',real=True)
m,alpha = sp.symbols('m,alpha',real=True,positive=True)

# 维度
d = int(input('输入维度: '))

if d < 2:
    raise ValueError('维度必须大等于2')

# 有心势
V = - 1/(d-2)*alpha/r**(d-2) if d >=3 else alpha*sp.log(r)

# 粒子Hamilton量
H = pr**2/(2*m) + ptheta**2/(2*m*r**2) + V

# 平衡轨道有dHdr1 = 0
dHdr1 = sp.diff(H,r).expand()
coeff_ptheta = dHdr1.coeff(ptheta,2)
coeff_others = dHdr1.coeff(ptheta,0)

# 平衡轨道的ptheta
# 理论解析解： sqrt(alpha*m*r**(-d+4))
_ptheta = sp.sqrt(- coeff_others / coeff_ptheta)
print(f'平衡轨道的ptheta:{_ptheta}')

dHdr = []
dHdr.append(dHdr1)

for i in range(1,5):
    dHdr.append(sp.diff(dHdr[-1],r))

print(f'H:{H}')
for i in range(0,5):
    print(f'd{i+1}Hdr{i+1}: {dHdr[i].subs(ptheta,_ptheta)}')
    # 理论解析解： dHdr = 0, d2Hdr2 = (4-d)*alpha/r**d
    # d2Hdr2 > 0代表稳定轨道(d=2, d=3)，否则是亚稳定(d=4)或者不稳定(d>4)
